global parity
global G
global H
global errTable
global n, 
global k
global n_k;

// (15, 7) double-error correcting code, generation polynomial X^8 + X^4 + X^2 + X + 1
n = 15;
k = 7;
n_k = n - k;
pol = [1 1 1 0 1 0 0 0 1]; 
[H, G] = cyclgen(n, pol);

errTable = syndtable(H);

H = H';
errBitNumber = 2;

function c = linearEncode(m)
  global G;
  c = modulo(m*G,2);
endfunction

function m = linearDecode(c)
  global H, errTable, k;
  len = size(c,1)
  m = zeros(len , k)
  for j = 1:len
    row = c(j, :);
    s = modulo(row * H,2);
    err = errTable(bi2de(fliplr(s))+1, :);
    row = modulo(err + row,2)
    m(j, :) = row(1:k);
  end
endfunction

EbN0_db = 0:11;
Eb = 1;       //signal energy per bit
EbN0 = 10^(EbN0_db./10);
N0 = Eb./EbN0;         //noise energy
sigma = sqrt(N0./2);

totalBitNum = k*1D7;
blockNum = 10000;
errBitNum = zeros(EbN0_db);
msgBer_157cyclic = zeros(EbN0_db);


for j = 1:length(EbN0_db)
  thisSigma = sigma(j);
  for i = 1:totalBitNum/(blockNum*k)
    msg = round( rand(blockNum, k, 'uniform')); 
    code = linearEncode(msg); //no coding
    code = 2*code - 1;
    noise = thisSigma*rand(blockNum, n, 'normal');
    rcv_code = sign(code + noise);
    rcv_code = (rcv_code + 1)/2;
    rcv_msg = linearDecode(rcv_code);
    errBitNum(j) = errBitNum(j) + length(find(rcv_msg ~= msg));
  end
end

msgBer_157cyclic = errBitNum./totalBitNum
